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Summary 


Nowadays models for self-organisation are being used in systems with a great 
degree of complexity and across disciplines. We show that the used Turing 
model is sensitive to parameters, type of domain growth, but also to the 
precision of model formulation itself. Hence it is necessary to revise Tur- 
ing’s model for self-organisation. For this purpose we consider derivation of 
evolution equations within non-equilibrium thermodynamic to identify phys- 
ically relevant formulations. Only then we subject these models to a detailed 
mathematical analysis. We offer possible extensions of the concept of self- 
organisation to more general situations and discuss its physical interpreta- 
tion. 


The essence and importance of these ideas is illustrated on the reaction- 
diffusion-advection system, where we indicate that such a system should be 
preferred from both physical and mathematical viewpoint. Further we point 
to the indispensable role of physical viewpoint during relevant model for- 
mulations. Using the non-equilibrium thermodynamic framework physically 
consistent extensions of Turing model are revealed as well as functional con- 
straints for present parameters. 


Souhrn 


V současné době se modely pro prostorové uspořádání hojně užívají i ve 
složitých systémech napříč disciplínami. Ukazujeme, že používaný Turingův 
model je citlivý na parametry, charakter růstu oblasti, ale i na samotnou přes- 
nost formulace modelu. Proto je nutné revidovat Turingův model pro vznik 
samovolného uspořádání. K tomuto účelu uvažujeme odvození evolučních 
rovnic v rámci nerovnovážné termodynamiky, abychom identifikovali fyzikálně 
relevantní formulace. Až poté podrobujeme tyto modely detailní matemat- 
ické analýze. Nabízíme též možná rozšíření konceptu sebeorganizace do obec- 
nějších situací a diskutujeme jejich fyzikální interpretaci. 


Podstatu a důležitost těchto myšlenek ilustrujeme na reakčně-difuzně- 
advekčním systému, kde naznačujeme, že takovýto systém by měl být prefer- 
ován jak z fyzikálního tak z matematického pohledu. Dále pak ukazujeme, 
jakou nezastupitelnou roli sehrává fyzikální pohled při tvorbě relevantních 
modelů, kdy pomocí konceptu nerovnovážné termodynamiky odhalujeme fy- 
zikálně konzistentí rozšíření Turingova modelu ale i funkční omezení pro 
vyskytující se parametry modelu. 
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1 Motivation 


Self-organisation in nature is widely recognised and is extensively modelled. 
Particularly, pattern formation due to chemical instability is believed to be 
of the essential importance in many non-equilibrium systems, ranging from 
developmental biology [19, 1] (e.g. morphogenesis [21], hair follicles [23] or 
fish pigmentation [31, 16]) to chemical reactions [3, 14] or in coupling of heat 
and mass transfer (the so-called Rayleigh-Bénard convection cells) [4]. For a 
review about reaction-diffusion models and a great variety of spatial patterns 
that they can exhibit see [5]. 

Systematic description of self-organisation in nature was initiated by Tur- 
ing [32] (from mathematical perspective) and also by Prigogine [26] (from 
non-eguilibrium thermodynamics perspective, NET). Turing showed that 
small local spatial fluctuations in an otherwise well-mixed system of auto- 
catalytic and inhibitory diffusing species (also known as morphogens) could 
become unstable due to diffusion and that an amplification of these fluctu- 
ations could lead to pattern development. Specifically, heterogeneous con- 
centrations of chemicals form a 'pre-pattern'. Subsequent differentiation of 
tissue/cell type is in response to whether or not concentration of one of these 
morphogens in the pre-pattern exceeds some threshold locally. Hence Tur- 
ing realised that this so-called diffusion-driven instability (DDI, or Turing 
instability, TI) can be considered as a symmetry breaking mechanism. Inter- 
estingly, Prigogine defines dissipative structures as a self-organisation that 
results from a connection between stability and dissipation [26] and where a 
pattern is driven by the distance from an equilibrium or nonlinearity in the 
system. Turing instability is, therefore, a special case of dissipation structures 
where particular choices of dissipation and system complexity are considered. 

For completeness we should mention that not all processes leading to pat- 
tern formation are of Turing type. To demonstrate this we take a well-studied 
example from developmental biology, the mechanism of spatial organisation 
in the Drosophila embryo. In this system a morphogen called bicoid is pro- 
duced at the future anterior pole of the embryo and diffuses from this source 
to form a concentration gradient across the embryo. This concentration gra- 
dient of bicoid ultimately serves as a positional guide, instructing nuclei of 
their position relative to the anterior pole. In the Drosophila embryo no true 
symmetry ever exists to be broken, as bicoid localisation to the future an- 


terior pole is achieved during the formation of the egg by maternal inputs 
[33, 7]. Thus, this morphogen gradient system leads to positional information, 
but the spatial pattern is not self-orchestrated. 

On the other hand, Turing instability represents a fundamentally different 
mode of pattern formation, and one which is capable of breaking symmetry 
without pre-existing positional information. This mechanism has driven nu- 
merous experimental studies even in the context of developmental systems 
(e.g [22, 8, 30, 24|) which suggest that Turing-like morphogen interactions 
and patterns can occur in such scenarios. However, a direct verification of 
Turing’s mechanism in biology at the level of molecular details has remained 
elusive. 

In this talk we shall first recall the classical work of Turing and identify 
some of its weaknesses. After that we will combine both the mathematical 
and physical viewpoints, draw ideas about plausibility of the used modelling 
approaches and also look for extensions there in order to assess a physically 
sound model for the emergence of self-organisation in nature. 


2 Diffusion-driven instability, Turing 


With a given model we use analytical tools to reveal when the self-organisation 
can be expected. ‘This is done via the so-called stability analysis. Its pre- 
cise meaning is context dependent and reflects what is meant by the self- 
organisation in a given system. Typically pattern formation is considered to 
correspond to a disruption of stability of a reference state due to a certain 
critical phenomenon taking place. Such a situation is also known as a bi- 
furcation and where the corresponding bifurcation parameter has a physical 
interpretation relevant for the considered context. In the case of standard 
Turing instability we shall discuss this matter in more detail below. 


2.1 Turing instability on static domain 


Diffusion-driven or Turing instability is defined for reaction-diffusion (RD) 
systems 


Ou = Dôu + f(u) 


with diagonal diffusion matrix D = diag( D1, D2) as a situation where a 
homogeneous steady state solution (HSS) u*, 0 = f(u*), is stable in the 


absence of diffusion and unstable once diffusion is introduced into the system. 
In this setting it can be equivalently rephrased as a requirement that the HSS 
is stable with respect to spatially homogeneous perturbations but unstable 
with respect to heterogeneous perturbations. As soon as we depart from such 
reaction-diffusion system with constant diffusion coefficients by considering a 
more complex transport operator, these two concepts may split up and also 
other admissible definitions of Turing instability may seem plausible. Hence 
we shall firstly discuss the motivation for the classical definition of Turing 
instability. For its discussion it is convenient to define a critical length and 
to relate it to TT. 

Let us solve a RD initial-boundary value problem in 1D with domain 
length L being a bifurcation parameter, i.e. to observe the effect of the (grow- 
ing) domain size on the stability and solution properties. Particularly, let us 
consider a one-dimensional RD model for a single species with, for example, 
Dirichlet boundary conditions 


Ou = D8 u + f(u), u(0,t) = 0 = u(L,t) 


where L > 0 is the size of the domain and with u* = 0 being the HSS, i.e. 
f(u*) = 0. To explore DDI we would like to solve the above system when 
accompanied by a small initial perturbation of the HSS 


0< u0) <1: 


Using the method of separation of variables and invoking completeness of 
eigenfunctions of Laplacian on L?(0, L) we have that the homogeneous solu- 
tion is of the form 


The nonhomogeneous solution is then obtained by expanding f(u) in terms 


of the orthonormal basis {sin (222) one As we are interested in small initial 
perturbations only (linear stability analysis), the expansion is at hand as 


fu) = f'(0)u(x). This gives that the time-dependent coefficients satisfy 


an 105 
For (linear) stability we require the solution to decay to the zero HSS or 
equivalently in this case all amplitudes of eigenfunctions have to vanish as 
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t — m. Thence we require that 


Dn? def 
L< EZO = Leori Vn. 
Thence there is a critical lengthscale Leri below which spatial self-organisation 
cannot occur as no eigenfunction fits in such a domain. 

These ideas extend to higher dimensions for Dirichlet boundary conditions 
as from the spectral properties of minus Laplacian for bounded domains it 
follows (from the method of Rayleigh quotient and minmax principles) that 
D0 C Q => An(Q1) > An(Q2) Vn with A, being the countable eigenvalues. 
The same conclusion does not hold for Neumann boundary conditions in 
general but holds for scaled domains, i.e. Q1 C Q2, Qo = aQy > An(Q1) > 
An(Q2), Vin. 


2.2 Why considered as a model for self-organisation? 


As we are about to study the emergence of self-organisation in more general 
systems, we first need to extend the concept of Turing instability to systems 
outside of RD. Two immediate possibilities from the above characterisation 


of TI are: 


S-O the emergence of self-organisation. We believe/propose that this situa- 
tion can be characterised by: one, the existence of a critical domain size 
below which system is stable to small perturbations and above which in- 
stability may occur; two, the existence of only a finite number of unstable 
modes as otherwise continuum description would break down (and no 
structure in spatial organisation would appear except salt-and-pepper 
pattern). This S-O instability can be also termed as a domain-size driven 
in which case the instability induces symmetry breaking (in a system 
that is continuously subjected to perturbations). 


TDI Transport-driven-instability requiring stability (of the HSS) without 
transport and instability once transport is considered. Naturally, the 
condition on the existence of finite numbers of unstable modes is present 
as well. Hence this TDI is equivalent to heterogeneity amplifying insta- 
bility, i.e. the system is stable until a small heterogeneity perturbs it, and 
thus the instability amplifies symmetry breaking rather than inducing 
it. 


Particularly, if one would like to study the transport-driven instability, one 
should adopt the above TDI concept of TI extension requiring the stability of 
the HSS without transport. Finally, if one addresses the spontaneous emer- 
gence of the spatial self-organisation in a reaction-diffusion-advection (RDA) 
system, one could adopt the above S-O extension which induces symmetry 
breaking (and hence self-organisation) while requiring the stability of the 
HSS only if zero is an eigenvalue of the transport operator with given BCs. 

Is it reasonable to promote diffusion-driven instability over other po- 
tential causes of instabilities when analysing the possibility of spatial self- 
organisation in a system? Why would one enforce the strict condition for 
stability of the HSS without diffusion in a system when this is not relevant 
to a real system where diffusion either is or is not present? In the case of 
RD system, answer to the first question is that the emergence of spatial self- 
organisation is related exactly to DDI as can be seen from the above S-O 
characterisation of Turing instability. Particularly, the requirement of the ex- 
istence of a critical length is related to the requirement of stability of the 
HSS in the following way. In the case of 


Neumann BCs the stability of a homogeneous steady state is required. The 
reason is that the existence of critical length enforces decay of perturba- 
tions for small enough lengths and hence, as orthonormal basis (ONB) 
contains a constant eigenfunction, stability of the HSS with respect. to 
spatially homogeneous perturbations is required. 


Dirichlet BCs the stability of the HSS is not implied from the stability of 
the system to perturbations as a constant is not part of ONB (nonzero 
spatially homogeneous perturbation does not satisfy Dirichlet BC). For 
example consider a 1D domain Q = (0, L], then eigenvalues of minus 
Laplacian are "F for a natural number n. As zero eigenvalue is not 
permissible due to the BCs, the zero homogeneous state is automatically 
stable in a linearised system for small enough L as eigenvalues of the 
linearised kinetics will be dominated by —"F [12] and hence without any 
condition on the stability of the HSS without diffusion. 


Therefore the stability of the HSS without diffusion in Turing instability (TT) 
is an additional requirement only for Dirichlet BC. 

Note that these considerations refer to a more general question: should 
Turing instability be a method/concept/model for an emergence of spatial 
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self-organisation or a particular type of (transport) driven spatial organisa- 
tion? Also note that the requirement of the existence of a critical length or 
that there is only a finite number of unstable modes might not be always 
required although they seem to be reasonable requirements for studying self- 
organisation in nature, especially in biology. In RD systems, these scenarios 
coincided as discussed above but in other situations we should choose an 
appropriate extension following a given aim. 


3 Is Turing model a plausible model for self-organisation? 


We shall now raise some of the known issues of Turing model and discuss 
whether it stands as a plausible model for self-organisation. 


3.1 Dı £ Də and binding to substrate 


One of the well-known problems of Turing’s approach is the requirement for 
unequal diffusion coefficients where for typical parameter values of reaction 
kinetics the diffusion coefficients of a putative morphogen pair need to dif- 
fer by an order of magnitude. This condition is, however, in contrast with 
Einstein-Smoluchovski relation for an estimate of diffusion coefficient as it 
would require the putative morphogens to significantly differ in size which is 
not observed. 

We recently showed that if at least one of the morphogens interacts with 
(binds to) underlying substrate (like the extracellular matrix), this issue can 
be resolved [17]. Consider linearised kinetics of two morphogens where one 
binds to a substrate 


Ou = D, Aut (fu — klu + fou + kw 
Ov = DAV guu + gyv 


Ow = ku — kw 


and Neumann BC for u, v. Note that k- > 0 (self-inhibitory) in order to 
have a chance for DDI [12]. 

Does fast binding change DDI conditions? Not in the quasi steady state 
limit as all the binding contributions dissapear. Instead we shall try to iden- 
tify a higher order approximation of this quasi steady state approximation. 


11 


Solving the last equation for bounded morphogen w yields (Laplace method): 
t t 
w(t)=e f kye- tu(r)dr = k, | eC —)u(r)dr = 
0 0 


t 
N F e 5- CDy(r)dr = 
t—e 
== ek- ek. = nek- | ko 
sky oz + amy a 


where the expansion of the solution u(t) = u(t) + (T — t)O;u(t) +... was 
used. Hence (ek_ < 1) 


w(t)» = | ut) dni) o (zz) 


Therefore the higher order quasi steady state approximation results in a 
rescaling of the relation for u 


(1 + =] Ou = D Aut fuu + fov 


>1 


Ov = D Av + guu + guv. 


Consequences on DDI conditions are easily assessed (the rescaling denoted 


by a prime (V affects D, and fu, fy): 


trJ’ < 0,.detJ > 0, Dugo + Do fa > 2V Da Dy det J > 0 


where D, > D, is no longer a contradiction. 


3.2 Reductionism 


As the above note about the prominent effect of considering binding of mor- 
phogens to a substrate indicates, the level of detail included in a model for 
self-organisation may significantly change the prediction of pattern forma- 
tion. 

The problem of reductionism was addressed in [12] pointing to the fact 
that a complete knowledge of morphogen interaction network is necessary 
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for plausible predictions of system’s behaviour. Particularly we showed that 
model predictions can be exactly opposite from reality. In one example the 
full model (say reality) does not yield patter but a reduced model, where one 
of the interacting morphogens is neglected (is undetected), does. Similarly, 
another example is provided but where the situation is reversed. 


3.3 Growing domain vs L as bifurcation parameter 


In DDI mechanism for the emergence of self-organisation the domain size 
is treated as a bifurcation parameter, 1.e. without the explicit time depen- 
dence. In 2010 Madzvamuse, Gaffney and Maini showed that tracking the 
actual slow growth has a qualitative impact on the results of analysis [18], 
e.g. a new type of reaction kinetics is allowed to give rise to DDI (although 
the parameter space is somewhat limited), different types of growth yield 
different Turing spaces (regions in parameter space where the necessary DDI 
conditions are met). Hence these findings are highlighting a qualitative differ- 
ence between the two approaches although the actual effect of domain growth 
on DDI conditions was found to be only via a correction to the classical DDI 
conditions of a small order. 

Based on this study we recently showed [13] that the actual picture of con- 
ditions for the emergence of a pattern on growing domains is more complex 
than previously thought and that the classical approach with domain size 
being the bifurcation parameter is a plausible approximation only for very 
slow growth. Additionally, for faster growth we show that prediction becomes 
even more out of reach. Particularly due to the intrinsic growth, modes can 
first substantially transiently grow prior to decay, but this behaviour is highly 
dependent on initial conditions or noise thence making predictions impossi- 
ble. Further, for some growth rates we observed that all modes with a high 
enough wavenumber transiently grow yielding a breakdown of the continuum 
description itself. 

For illustration consider a RD system on a growing domain which can be 
simplified to the following dimensionless formulation 


ES Cre sap Poe +4F(u) (1) 


in normalised Lagrange coordinates! and with the following definitions of 


Ir = x(E,t) where y(€,0) = € and č € [0,1] 
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parameters: D = diag(1,d), d = D,/Du, y = L2/Dy. The advection term, 
V.(au), representing growth is transformed into an effective time-dependent 
diffusion coefficient”, nD, and an additional linear term h(t)u. This ex- 
plicit dependence of the RD equation on time as a consequence of intrinsic 
domain growth is making the stability analysis difficult. However, with the 
assumption of a slow uniform linear growth one can approximate the qual- 
itative evolution of a perturbation about a reference state u, (a spatially 
homogeneous but time-dependent solution to the RD problem) yielding an 
analogue to DDI conditions for pattern formation [18]. 

A point in parametric space is a member of the so-called Turing space 
iff the corresponding reference state is stable without diffusion and unstable 
once diffusion is present for the reasons discussed above. Hence a point in 
the parametric space is not in Turing space iff either the reference state itself 
is unstable or that diffusion failed to destabilise the stable reference state. 
In the classical treatment one can access this information via evaluation of 
DDI conditions for such a value of parameters. But when DDI conditions 
are evolving, as is the case with the intrinsically growing domain, Turing 
space has to be more carefully assessed. As a result, when time-evolving DDI 
conditions are plotted (as reported in [18] e.g. for linear growth) one needs 
to track the reason why a given point changed its property with respect to 
DDI conditions. In fact, if the plotted DDI region is moving into a region 
of parameter where the reference state itself was unstable (being the reason 
why Turing space was not here in preceding times) this region cannot be 
considered as DDI leading to self-organisation (e.g. due to the lack of exis- 
tence of critical length) and, as a result, should not be included in the Turing 
space. Hence, one needs to track the evolution of stability of the reference 
state in addition to evolution of intersection of DDI conditions in order to 
correctly capture the parametric space leading to the emergence of spatial 
self-organisation on growing domains. 

Particularly, Madzvamuse et al [18] showed that DDI conditions evolve as 


2x(€,t) = plt)E = €exp( Ji h(r)dr) 
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follows 


=Y tr Tasta) + 2h(t 
—h(ti)y tr Jan) + 7° det n (t1) 


=q [dfu+ gu] + h(ti)(1 +d) <0 
[A(t ) (1 + d) m y(dfu + 20 = 4dlý det Ju, (tı) m yh(t) tr Jae) > 0 

U(t) 

where Ju) denotes the Jacobian of linearised kinetics F(u) evaluated at 


u,(t;) and fu, fo, Ju ,Gv represent its elements. Consider t € [0,7]. A 
bifurcation due to a DDI occurs, for the first time, at t = T > 0 on the set 


Use| nudí Uue \ [Sumo]. 
te[0,7) n te(0,T) —y— 
at time T —v assuming not 
stability w.r.t. homogeneous but not unstable i ra 
perturbations always at any earlier time t>0 at I= 


Turing instability will have occurred in the ‘Turing space below by time T 


sonzojví U || Us@)a(ao- U Aw 
=-€V—V—ae—~—_—__—" s€(0,T] te [0,5] te|0,s] 
DDI occurs at t=0 


Turing instability for first time at t=s 
neglecting subset where in fact DDI occurs at t=0 


Note that these ideas were illustrated for two types of reaction kinetics 
where the complexity of conditions yielding self-organisation is apparent [13]. 


3.4 NET - constantness of D? 


Physical perspective and insight is a key in mathematical modelling. As we 
shall illustrate in this and the following subsections, a detailed mathematical 
analysis of a given problem (in this case the Turing’s RD problem with con- 
stant diffusion coefficient) may yet provide misleading conclusions due to the 
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weak correspondence of the model to reality. In all the three above examples 
point to the sensitivity of the results to the model formulation and hence the 
formulation of a model for self-organisation should be careful. 

We shall show here that the typical assumption of constantness of diffusion 
coefficient in Fick’s law may be in many situations false by yielding a thermo- 
dynamically inconsistent model. We use the non-equilibrium thermodynamic 
framework for formulating physically sound models. A typical and most used 
example of such framework is the classical irreversible thermodynamics (CIT) 
which has been shown to be a very powerful tool for modelling transport, re- 
action kinetics or viscous effects. The theory has been described [20] and 
thoroughly summarized |6]. It has been applied to engineering problems, in- 
cluding coupled heat and mass transport |10]. A fundamental assumption of 
the theory are the linear force-flux relations 


Ji = a Lis Xj, (2) 
J 


where J; and X; are the thermodynamic fluxes and forces, respectively. Phe- 
nomenological coefficients L;; form a matrix of phenomenological coefficients 
L and in general they depend on the local thermodynamic state of the system, 
expressed in terms of the state variables. 

The success of CIT is closely related to the accessibility of experimental 
assessment of the closure (constitutive) relations. Nevertheless, for n coupled 
processes n? phenomenological parameters need to be identified together with 
their dependence on state variables (and boundary conditions). Not all are, 
however, easily accessible and hence any a-priori relations among the co- 
efficients mean a significant reduction of the necessary experimental effort 
to assess all the model parameters. To this end the second law of thermo- 
dynamics, which is expressed as non-negativity of local entropy production 
os = >), J;X; > 0, imposes constraints on L, particularly L being positive 
semi-definite, and the Onsager reciprocal relations (ORR) [27, 28] reduce the 
number of unknown parameters to n?/2 + n/2. 

Recently we showed |15] that there are further constraints on phenomeno- 
logical coefficients in addition to ORR within CIT where this claim was 
motivated by our observation about correlation of transport coefficients in 
water and proton transport across Nafion membranes at low applied elec- 
tric potentials [2]. We hypothesised that in addition to Onsager’s relations 
the coefficients have to share the same functional dependence on the local 
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thermodynamic state of the system and proved this statement for a certain 
class of nonequilibrium thermodynamic models using the maximum principle. 
Thence mathematical models and experimental data should be checked for 
consistency in the above sense and we discussed this issue on some examples. 
Particularly, we showed that (by studying a thermodiffusion model which in 
an isothermal case reduces to Fick’s diffusion law) that the widely used as- 
sumption of constant diffusion coefficient is thermodynamically inconsistent 
once the considered system would have a nonzero thermodiffusion coupling 
in the nonisothermal case, i.e. the presence of temperature gradients would 
cause a diffusion flux to appear as in the classical Rayleigh-Bénard convec- 
tion. On the other hand, if we let the diffusion coefficient to scale as D œ c~! 
with c being the concentration we obtain a thermodynamically consistent 
model. 


3.5 NET (theory of mixtures) - extension of RD to RDA 


If we subject the Fick’s law of diffusion to scrutiny, particularly if we derive 
the evolution equations for reacting mixtures from first principles and use the 
CIT or EIT (extended irreversible thermodynamic |9]) constitutive theory, 
we observe that advection can appear as a pure consequence of chemical reac- 
tions among constituents. Particularly, from a careful formulation of mixture 
theory of fluids [29] it follows that chemical kinetics is not only driven by 
chemical affinity A, but by diffusion flux ug as well. We have shown that the 
chemical reaction rate E, is governed by 


: rr 1 a 
Či = B Lrs (a + 2 nb) ; 


where Vsa denotes stoichiometric coefficient of species a in s-th reaction and 
M, stands for molar mass. The effect of affinity on chemical reaction rate 
is typically modelled by Law of mass action. Crucially, the identified term 
reveals a new effect: diffusion can drive chemical reactions but also vice-versa, 
chemical reactions alone can cause diffusion flux. 

Substitution of this relation together with Fick’s law of diffusion into the 
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balance of mass yields 


Oia = V(DaVea) + X VraMabr = DaV?ea + X Ka(Veg)? + R(ei,..., en) 
r B 


= DgVea + X Ka(Veg) (Veg) + R(c,...,(n), 
B M 


where the kinetics R(c1,..., Cn) is determined from law of mass action and 
interaction network |11]. Note that the new term may be formally regarded as 
an advection. Moreover, advection with velocity Vy may be present regardless 
of the presence of Fickian diffusion and thus should be taken into account. 
In summary, this suggests that the effect of (small) advection on RD models 


should be studied. 


3.6 Self-organisation in RDA systems 


Let us consider a RDA system and for simplicity let us restrict the analysis 
to a 1D domain G = [0, L]. Particularly, we have 


r(z)a\ — (Da L(a) f (a,b) 
a ts i) s (a eo) 7) a,b) 
where La = 9, (p(x)0,a) + q(x)a is an operator discussed below, Da, Dp are 
two constants and a,b are concentrations of two morphogens. 
First, note that stability analysis of a stationary solution that is spatially 


dependent is beyond the scope of this work. Hence, we restrict ourself to the 
standard case when the existence of a homogeneous stationary solution a*, b* 


M Pa f(a*, bř) 
the fraction g(x)/r(x) is a constant and hence the q(x) term can be absorbed 
into the reaction kinetic terms. As a result we can assume WLOG q(x) = 0 


is required. Thus 


in the subsequent analysis, although such a shift of linearised kinetics can 
have a profound yet straightforward effect on the character of stability of the 
HSS and also on the emergence of self-organisation. 

Hence we consider the above equations subjected to non-homogeneous 
Robin, a1y(0) + agy’(0) = asy, Piy(L) + Poy'(L) = Biy*, or periodic, 
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y(0) = y(L), y'(0) = y'(L), boundary conditions that are compatible with 
the HSS y*. The spatial operator £ = 9, (p(x)0,) with these boundary condi- 
tions is a Sturm-Liouville operator. It is a self-adjoint operator when defined 
on L?(0, L) and the set of its eigenfunctions satisfying the above bound- 
ary conditions form a complete basis set of those sguare-integrable functions 
which satisfy the given boundary conditions as can be shown with the help of 
Rayleigh quotient. Note that for the periodic BCs they also form a complete 
set of eigenfunctions, only in this case there might not be a unigue linearly 
independent eigenfunction corresponding to each eigenvalue and would need 
to use Gramm-Schmidt process first. 
Let n(x) solve 0, (p(x)O,n(x)) = —Anr(z)yn(r) and let 


OROUA 


with (3) = |, AnYn(x) where we take advantage of the fact that yn(£) 


1 
form an ON basis of L? (G) and A, = i) 


b* 


So r(2)m(2) (a, + Coe — J(a*, PAn ) -0 


n 


Linearising about ) yields 


where J(a*, b*) denotes the jacobian of reaction kinetics evaluated at the HSS 
(e.g. Jig = Ohfla=a*,b=b+). AS Yn form ON basis we have that 


Da 0 
An=- On (75 p) 1) A 


and hence the so-called dispersion relation is obtained 


D0 7 
di (outdo (Pe 8) 3) 0 


Its roots on correspond to eigenvalues of the linearised system governing the 
linear behaviour and consequently determine the linear stability of the HSS. 

The difficulty, of course, lies in the eigenvalue problem of the S-L operator 
L. To illustrate the above ideas we shall discuss in detail a situation of a RDA 


19 


problem and study what implication it has on TDI and S-O in relation to 
Turing instability. We shall consider both possible concepts and extensions 
of Turing instability proposed above and study systems with advection. The 
main aim of this section is to estimate the consequences of negligence of 
(small) advection in Turing instability studies. 

We consider a system of two interacting species a, b without differential 
transport, i.e. where both species are diffusing and advecting with the same 
magnitude (Da = Dy and advection is the same) 


a D0?,a + VO,a f(a, 5) 

à (5) = Coe a) n Cae (3) 
where functions f,g describe reaction kinetics. The above S-L theory can 
be applied as the choice r(x) = p(x) = e“"/P with D = D, yields a RDA 
system with a constant advection V. 

Three types of boundary conditions are relevant for RDA problems. Par- 


ticularly with (a*, b*) being the HSS about which linear stability is of interest 
we explore (analogous BCs are required for species b) 


Dirichlet BC a(x) = a* at x € 10, L), 

Fixed-flux BC D,0,a(x) + Va(x) = Va* at x € {0, L}, 

Periodic BC a(0) = a(L) and 0,a|,=9 = O,a|:=L, 

Danckwerts BC D,0,a(x) + Va(x)|,-1 = Va* and 0,a(x)|,—9 = 0, 

so that the HSS is a plausible solution of the problem at hand. However, for 
the sake of illustration, we shall consider only the fixed-flux BC for species 
a, b (other BCs can be analogically analysed), i.e. La = Va(x)+D0,a(x) = 
Va* and £% = Vb* at x € {0, LY, and a small perturbation of the HSS as 
an initial condition, i.e. a(t = 0,x) = a* + ao(x), b(t = 0, x) = bř + bo(z), 
lao| < a*, [bo] X bř. Compatibility of the boundary and initial conditions 


requires Vao(x) + DO,ao(x) = 0 at the boundary and similarly for species b. 
Linearisation about the HSS (a*, b*) with (a,b) = (a,b) — (a*, b*) yields 


TERORO : 


under the assumption of small perturbations (aiming at linear stability anal- 
ysis and hence dropping the higher order terms), |a| < a*,|b| « bř, and 
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with the matrix of linearised kinetics, J, evaluated at the HSS (a*, b*). The 
nonhomogeneous Robin boundary condition is transformed into a zero-flux 
boundary condition (L°a, C*“b) = 0 at both boundaries of [0, L]. The initial 
conditions are a(t = 0, £) = ao(x), b(t = 0, x) = bo(x) and are automatically 
compatible with the zero-flux boundary conditions. 

Now we are ready to employ the ONB of L?(0, L) that consists of eigen- 
functions yn of L, i.e. Cyn = An Yn subjected to homogeneous BCs L°°y,, = 0, 
as all the terms appearing in the bulk equation satisfy zero-flux BCs. Hence, 
stability of the HSS (a*,b*) can be assessed from the decays of the am- 
plitudes of the eigenmodes in the problem (4). With a = >, An(t)yn(2), 


b = >, Br(t)yn(x) we have 


` la ee (bd) 6] Yn = 0. 


n 


Thus the HSS is linearly stable if the roots ox of the dispersion relation 
det (—o+I — (A,I— J)) = 0 (5) 


have negative real parts and is unstable if at least one of the roots has a 
positive real part. 
Except for fine parameter tuning (that would be hard to justify due to the 
; 0 l 
required robustness) we have J = UT ) U and hence the dispersion 
u- 


relation (5) can be rewritten as 


det fur (-01 + (= ~An ° )) u| = 0 
0 U — An 


with roots o+ = u+ — An. Thus we may conclude (it can be easily shown 
that the conclusion does hold even when J is not diagonalisable via Jordan 


blocks) that the HSS (a*, b*) is linearly stable iff 


R(An — u+) > 0 and RIA, — w_) > 0, Vn E No (6) 


where A, are eigenvalues of the spatial operator £ subjected to BCs £L°°u = 0 


at x € {0, L} and u+ are two eigenvalues of linearised kinetics matrix J. 
To conclude the stability analysis it suffices to identify eigenvalues of the 
spatial operator £. Namely, 


= — (D97, + Vô) Yn = Ann 
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subjected to zero-flux BCs 
Ln = (Doet Vym 0 atr Ee 0L} 


The spatial operator £ is a non-self-adjoint operator but can be transformed 
into a self-adjoint operator [25] via transformation m = eD On with appro- 
priately modified BCs (note that the eigenfunctions €, are orthogonal with 
respect to the weighted inner product (£n, 6x) = Te €,(x)&,(v)eD"dax = 0 for 
n Æ k). We shall use, however, the above outlined S-L theory to illustrate its 
generality and applicability. With the choice of p(x) = r(x) = e"“/P, finding 
eigenfunctions requires solving 


6 a = Xn 


with constraints Vy, + Dôn = 0 at x € 40, L}. One finds that such 
eigenvalues are 


nr? 1V? 
0. No. 
a ID AE 


Hence we may summarise this subsection that in the case of prescribed 
fixed-flux BCs there is no TDI as it requires stability of the HSS (and hence 


Rus < 0, Ru- < 0) but for (S-O) we may have an emergence of self- 
organisation for arbitrarily small magnitude of advection V where the fixed- 


M=DÍ 


flux BCs approach zero-flux BCs used in RD systems. This is in contrast 
with classical TI in RD systems but also in the case when we would adopt 
the definition of (S-O) for TI. The reason is that in the classical problem 
the stability of the HSS is required as zero-flux BCs allow for a constant 
(a HSS) solution (zero is an eigenvalue) whereas the fixed-flux BCs in the 
case of RDA do not permit the HSS to be part of the ONB (zero is not an 
eigenvalue) and hence stability of the HSS is not reguired for the emergence 
of self-organisation. Finally, note that the existence of critical length below 
which instability cannot occur together with the finite number of unstable 
modes are satisfied. 


4 Conclusion 


Let us summarise the main observations from this talk as 
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1. we do not have a plausible model for self-organisation (for applications 
in nature) 


2. nonequilibrium thermodynamics is a suitable framework to use to iden- 
tify model formulations 


3. only then (after a suitable model is identified) a detailed mathematical 
analysis is needed and can reveal conditions for the emergence of self- 
organisation. 


References 


[1] R.E. Baker, E.A. Gaffney, and P.K. Maini. Partial differential equations 
for self-organization in cellular and developmental biology. Nonlinearity, 
21(11):R251, 2008. 


J.B. Benziger, M.J. Cheah, V. Klika, and M. Pavelka. Interfacial con- 
straints on water and proton transport across nafion membranes. Journal 
of Polymer Science Part B: Polymer Physics, 53(22):1580—1589, 2015. 


P. Borckmans, G. Dewel, A. De Wit, E. Dulos, J. Boissonade, F. Gauffre, 
and P. De Kepper. Diffusive instabilities and chemical reactions. Int J 
Bifurc Chaos, 12(11):2307—2332, 2002. 


[2 


— 


[3 


— 


[4 


— 


S. Chandrasekhar. Hydrodynamic and hydromagnetic stability. Interna- 
tional Series of Monographs on Physics Oxford, England . Dover Publi- 
cations, dover ed edition, 1981. 


[5] M.C. Cross and P.C. Hohenberg. Pattern formation outside of equilib- 
rium. Rev Mod Phys, 65(3):851, 1993. 


[6] S.R De Groot and P. Mazur. Non-equilibrium thermodynamics. Courier 
Corporation, 2013. 


[7] S.F. Gilbert. Developmental Biology. Sinauer Associates Inc, Sunder- 
land, Massachusetts USA, 8th edition, 2006. 


[8] M.P. Harris, S. Williamson, J.F. Fallon, H. Meinhardt, and R.O. Prum. 
Molecular evidence for an activator-inhibitor mechanism in develop- 
ment of embryonic feather branching. Proc. Natl. Acad. Sci. USA, 
102(33):11734-11739, 2005. 


23 


[9] D. Jou, J. Casas-Vázauez, and G. Lebon. Extended irreversible thermo- 
dynamics. Springer, 1996. 


[10] S. Kjelstrup and D. Bedeaux. Non-Equilibrium Thermodynamics of Het- 
erogeneous Systems. Series on Advances in Statistical Mechanics. World 
Scientific, 2008. 


[11] V. Klika. Comparison of the effects of possible mechanical stimuli on 
the rate of biochemical reactions. J Phys Chem B, 114(32):10567—10572, 
2010. 


[12] V. Klika, R.E. Baker, D. Headon, and E.A. Gaffney. The influence of 
receptor-mediated interactions on reaction-diffusion mechanisms of cel- 
lular self-organisation. Bull Math Biol, 74(4):935—957, 2012. 


[13] V. Klika and E.A. Gaffney. On the differences in diffusion-driven in- 
stability on static and growing domains. Proc Roy Soc A, July 2016. 
submitted. 


[14] V. Klika and M. Grmela. Mechano-chemical coupling in belousov- 
zhabotinskii reactions. The Journal of Chemical Physics, 
140(12):124110, 2014. 


[15] V. Klika, M. Pavelka, and J.B. Benziger. Functional constraints on phe- 
nomenological coefficients. Physical Review E, July 2016. submitted. 


[16] S. Kondo and T. Miura. Reaction-diffusion model as a framework for 
understanding biological pattern formation. Science, 329(5999):1616— 
1620, 2010. 


[17] K. Korvasová, E.A. Gaffney, P.K. Maini, M.A. Ferreira, and V. Klika. 
Investigating the turing conditions for diffusion-driven instability in the 
presence of a binding immobile substrate. Journal of theoretical biology, 
367:286—295, 2015. 


[18] A. Madzvamuse, E.A. Gaffney, and P.K. Maini. Stability analysis of non- 
autonomous reaction-diffusion systems: the effects of growing domains. 
Journal of mathematical biology, 61(1):133-164, 2010. 


[19] P.K. Maini, K.J. Painter, and H.N.P. Chau. Spatial pattern forma- 
tion in chemical and biological systems. J Chem Soc, Faraday Trans, 
93(20):3601-3610, 1997. 


24 


[20] J. Meixner and H.G. Reik. Thermodynamik der Irreversible Prozesse, 
in Handbuch der Physik, volume 3/II. Springer, Berlin Heidelberg New 
York, 1959. 


[21] R.J. Metzger, O.D. Klein, G.R. Martin, and M.A. Krasnow. The branch- 
ing programme of mouse lung development. Nature, 453(7196):745—750, 
2008. 


[22] T. Miura and K. Shiota. Extracellular matrix environment influences 
chondrogenic pattern formation in limb bud micromass culture: Experi- 
mental verification of theoretical models. Anat. Rec., 258:100—107, 2000. 


[23] C. Mou, B. Jackson, P. Schneider, P.A. Overbeek, and D.J. Headon. 
Generation of the primary hair follicle pattern. Proc Natl Acad Sci, 
103(24):9075-9080, 2006. 


[24] C. Mou, F. Pitel, D. Gourichon, F. Vignoles, A. Tzika, P. Tato, L. Yu, 
D. W. Burt, B. Bed'hom, M. Tixier-Boichard, K. J. Painter, and D. J. 
Headon. Cryptic patterning of avian skin confers a developmental facility 
for loss of neck feathering. PLoS Biol., 9(3):e1001028—, 03 2011. 


[25] K.N. Murty, L.V. Fausett, and D.W. Fausett. Transformation of a class 
of non-self-adjoint systems into self-adjoint hamiltonian systems. Jour- 
nal of mathematical analysis and applications, 174(1):1-21, 1993. 


[26] G. Nicolis and I. Prigogine. Self-organization in nonequilibrium systems. 


Wiley New York, 1977. 


[27] L. Onsager. Reciprocal relations in irreversible processes. I. Phys. Rev., 
37:405-426, Feb 1931. 


[28] L. Onsager. Reciprocal relations in irreversible processes. ii. Phys. Rev., 
38:2265—2279, Dec 1931. 


[29] M. Pavelka, F. Maršík, and V. Klika. Consistent theory of mixtures 
on different levels of description. International Journal of Engineering 
Science, 78:192—217, 2014. 


[30] L. Solnica-Krezel. Vertebrate development: Taming the nodal waves. 


Curr. Biol., 13:R7-R9, 2003. 


25 


[31] L Solnica-Krezel et al. Vertebrate development: taming the nodal waves. 


Current biology: CB, 13(1):R7, 2003. 


[32] A. Turing. The chemical basis of morphogenesis. Phil Trans R Soc Lond 
B, 237:37-72, 1952. 


[33] L. Wolpert. Principles of Development. Oxford University Press, 2nd 
edition, 2002. 


26 


Curriculum Vitae 
Ing. Vaclav Klika, Ph.D., born 12. 3. 1983 


Education 
October 2006 - October 2009: Czech Technical University in Prague, FNSPE, Special- 
ization: Mathematical engineering (Doctoral degree programme), Ph.D. thesis: Towards 
long-term prediction of tissue remodelling, supervisor: prof. Ing. F. Maršík, DrSc. 
September 2001 - June 2006: Czech Technical University in Prague, FNSPE, Special- 
ization: Mathematical modelling (Master degree programme), Masters thesis: Mathemat- 
ical model of bone remodelling, supervisor: prof. Ing. F. Maršík, DrSc. 


Work experience 
Occupation 
October 2009 - present: Assistant professor at department of Mathematics, FNSPE, CTU, 
Supervisor and co-supervisor of 2 undergraduate students and 4 graduate students 

January 2004 - December 2015: Researcher in the Institute of Thermomechanics (part- 
time), Academy of Sciences of the Czech Republic 
Projects 
EMBO short-term fellowship award for project “A systemic viewpoint on branching mor- 
phogenesis” placed at ETH Zürich, ASTF No: 123 — 2012 

Member of research team in grants: 
Modelling The Impact Of Collagen Anisotropy Within Cartilage (Platform grant of Math- 
ematical Institute, University of Oxford awarded to EA Gaffney); CENTEM project 
(CZ.1.05/2.1.00/03.0088) cofunded by the ERDF and through CENTEM PLUS (LO1402) 
from the Ministry of Education, Youth and Sports; Optimization of the chemical composi- 
tion, structural characteristics and mechanical properies of NiTi alloys for biomechanical 
applications (2009-12; GA106/09/1573); Material properties of veins and their remode- 
lation (2008-11; GA CR No. 106/08/0557); Materials and components for environment 
protection (2006-11; 1M06031); 


Stays abroad 
July - September 2015: A three-months stay at Newton’s Institute, Cambridge Univer- 
sity (research programme Coupling Geometric PDEs with Physics for Cell Morphology, 
Motility and Pattern Formation, 13.7 - 18.12. 2015) 

May-July 2014: A two-months stay in Centre of Mathematical Biology, Mathematical 
Institute, University of Oxford; Platform grant with EA Gaffney 

July-September 2012: A two-months stay at BSSE (D. Iber), ETH, Switzerland; project: 
Osteogenesis and branching morphogenesis as Turing instability on growing domains 

May 2011: A one-month stay at the OCCAM, Mathematical Institute, University of 
Oxford; continuation of initiated research; short-term visitor of OCCAM 

July - October 2010: A three-months stay at the OCCAM, Mathematical Institute, 
University of Oxford; Visiting Postdoctoral Research Assistantship award 

May - July 2010: A two-months stay in Research Group on Structural Mechanics and 
Materials Modeling (GEMM), I3A, University of Zaragoza; travel grant from Instituto 
Aragones de Ciencias de la Salud, research project PIPAMER10/015 

October 2005: One month stay in LIAB (prof Yahia) of Ecole Polytechnique, Montreal 


27 


Awards 

e the second place in Reinhart Heinrich Doctoral Thesis Award 2011 from European 
Society of Mathematical and Theoretical Biology 

e honourable mention in Votruba prize 2010 contest for the best thesis in theoretical 
physics from Doppler institute, Czech Republic 

e Prof. Valenta and Prof. Čihák award from Czech Society of Biomechanics 2007 
(master thesis) 


Other 
e active collaboration with EA Gaffney and PK Maini (Oxford), M Grmela (Mon- 
treal), Y. Kevrekidis and J. Benziger (Princeton), H Bougherara (Ryerson), D. 
Headon (Edinburgh), D Iber (ETH), JM Garcia-Aznar and M DoblarAT (Zaragoza), 
S Cotter (Manchester) 
e editorial and review work 
— editor of two volumes Theoretical Biomechanics and Biomechanics in Appli- 
cations by InTech, Vienna; August, November 2011, ISBN 978-953-307-312-5 
and 978-953-307-969-1 
— field Editor of General Physiology and Biophysics 
— reviewer of Journal of Chemical Physics, Journal of Mathematical Biology, 
Journal of Theoretical biology, Interface Focus, Entropy, International Journal 
of Engineering Sciences, Chaos, Physics Letters A, Mathematical Methods in 
the Applied Science, Biomath, Applied Bionics and Biomechanics, Journal of 
Thermal Science and Engineering Applications 
e member of a scientific committee for European Society of Biomechanics’ Congress 
2015 held in Prague, July 1-8 2015 
e contract with deGruyter publishing house, monograph "Non-eguilibrium Thermo- 
dynamics of Mixtures" (2016) 
e Scientometric data*; citations (WoK): 107 (67 without self-citation); H-index 6 
e invited lecture: 
— Mathematics of Pattern Formation 2016, Mathematical Research and Confer- 
ence Center, Bedlewo, Poland, September 12-17 2016 
— Keynote lecture, BIOMATH 2014, V Klika: Emergence of spatial organisation 
in real systems and its modelling, Sofia, Bulgaria, 22-27 June 2014 
— The European Conference on Mathematical and Theoretical Biology ECMTB 
2011, V. Klika, F. MarAaAnk: Tissue adaptation driven by chemo-mechanical 
coupling with application to bone in mini-symposium on Mathematical model- 
ing of biomechanical regulation in bone tissue, Krakow, Poland, June 28- July 
2 2011 
— V. Klika - Biochemical model of bone remodelling including mechano-chemical 
coupling. Banff International Research Station workshop, topic Mathematical 
Foundations in Mechanical Biology, September 2010 


3as of 6.7. 2016 


28 


10. 


11. 


12. 


13. 


List of publications 


. Pavelka, M., Klika, V. Esen, O., Grmela, M. (2016). A hierarchy of Poisson brackets 


in non-equilibrium thermodynamics. Physica D: Nonlinear Phenomena. Accepted 


. Benziger, J. B., Cheah, M. J., Klika, V., Pavelka, M. (2015). Interfacial constraints 


on water and proton transport across nafion membranes. Journal of Polymer Science 
Part B: Polymer Physics, 53(22), 1580-1589. 


. Grmela, M., Klika, V., Pavelka, M. (2015). Reductions and extensions in mesoscopic 


dynamics. Physical Review E, 92(3), 032111. 


. Korvasová, K., Gaffney, E. A., Maini, P. K., Ferreira, M. A., Klika, V. (2015). 


Investigating the Turing conditions for diffusion-driven instability in the presence 
of a binding immobile substrate. Journal of theoretical biology, 367, 286-295. 


. Pavelka, M., Klika, V., Vágner, P., Maršík, F. (2015). Generalization of exergy 


analysis. Applied Energy, 137, 158-172. 


. Avval, P. T., Samiezadeh, S., Klika, V., Bougherara, H. (2015). Investigating stress 


shielding spanned by biomimetic polymer-composite vs. metallic hip stem: A compu- 
tational study using mechano-biochemical model. journal of the mechanical behavior 
of biomedical materials, 41, 56-67. 


. Cotter, S. L., Klika, V., Kimpton, L., Collins, S., Heazell, A. E. (2014). A stochas- 


tic model for early placental development. Journal of The Royal Society Interface, 
11(97), 20140149. 


. Pavelka, M., Klika, V., Grmela, M. (2014). Time reversal in nonequilibrium ther- 


modynamics. Physical Review E, 90(6), 062131. 


. Klika, V., Pérez, M. A., García-Aznar, J. M., Maršík, F., Doblaré, M. (2014). A 


coupled mechano-biochemical model for bone adaptation. Journal of mathematical 
biology, 69(6-7), 1383-1429. 


Avval, P. T., Klika, V., Bougherara, H. (2014). Predicting bone remodeling in re- 
sponse to total hip arthroplasty: computational study using mechanobiochemical 
model. Journal of biomechanical engineering, 136(5), 051002. 


Pavelka, M., Maršík, F., Klika, V. (2014). Consistent theory of mixtures on different 
levels of description. International Journal of Engineering Science, 78, 192-217. 


Klika, V., Grmela, M. (2014). Mechano-chemical coupling in Belousov-Zhabotinskii 
reactions. The Journal of chemical physics, 140(12), 124110. 


Klika, V. (2014). A guide through available mixture theories for applications. Critical 
Reviews in Solid State and Materials Sciences, 39(2), 154-174. 


29 


14. 


15. 


16. 


17. 


18. 


19. 


20. 


Klika, V., Grmela, M. (2013). Coupling between chemical kinetics and mechanics 
that is both nonlinear and compatible with thermodynamics. Physical Review E, 
87(1), 012141. 


Klika, V., Baker, R. E., Headon, D., Gaffney, E. A. (2012). The influence of receptor- 
mediated interactions on reaction-diffusion mechanisms of cellular self-organisation. 
Bulletin of mathematical biology, 74(4), 935-957. 


Bougherara, H., Klika, V., Maršík, F., Mařík, I. A., Yahia, L. H. (2010). New pre- 
dictive model for monitoring bone remodeling. Journal of Biomedical Materials Re- 
search Part A, 95(1), 9-24. 


Klika, V., Maršík, F. (2010). A thermodynamic model of bone remodelling: the 
influence of dynamic loading together with biochemical control. J. Musculoskeletal 
Neuronal Interact, 10(3), 220-230. 


Klika, V. (2010). Comparison of the effects of possible mechanical stimuli on the 
rate of biochemical reactions. The Journal of Physical Chemistry B, 114(32), 10567- 
10572. 


Maršík, F., Klika, V., Chlup, H. (2010). Remodelling of living bone induced by 
dynamic loading and drug delivery-Numerical modelling and clinical treatment. 
Mathematics and Computers in Simulation, 80(6), 1278-1288. 


Klika, V., Maršík, F. (2009). Coupling effect between mechanical loading and chem- 
ical reactions. The Journal of Physical Chemistry B, 113(44), 14689-14697. 


Chapters in books 


1. 


Klika V, Maršík F, August 2011. Biomechanics, Theory. IN-TECH, Vienna, Ch. Fea- 
sible Simulation of Diseases Related to Bone Remodelling and of Their Treatment, 
ISBN 978-953-307-312-5. 


. Klika V, Maršík F, Mařík I, February 2010. Dynamic Modelling. IN-TECH, Vienna, 


Ch. Influencing the Effect of Treatment of Disease Related to Bone Remodelling by 
Dynamic Loading, ISBN 978-953-7619-68-8. 


. Bougherara H, Klika V, Maršík F, Mařík I, Yahia L, 2009. Damage and Fracture 


Mechanics. Springer Netherlands, Ch. A Novel Approach for Bone Remodeling After 
Prosthetic Implantation, pp. 553-565, ISBN 978-90-481-2668-2. 


30 


